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We derive a model describing spatio-temporal organization of an array of microtubules interacting 
via molecular motors. Starting from a stochastic model of inelastic polar rods with a generic 
anisotropic interaction kernel we obtain a set of equations for the local rods concentration and 
orientation. At large enough mean density of rods and concentration of motors, the model describes 
an orientational instability. We demonstrate that the orientational instability leads to the formation 
of vortices and (for large density and/or kernel anisotropy) asters seen in recent experiments. We 
derive the specific form of the interaction kernel from the detailed analysis of microscopic interaction 
of two filaments mediated by a moving molecular motor, and extend our results to include variable 
motor density and motor attachment to the substrate. 
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I. INTRODUCTION 

One of the most important functions of molecular mo- 
tors is to organize a network of long filaments (micro- 
tubules) during cell division to form cytosceletons of 
daughter cells [lj . In order to better understand the de- 
tails of this complex self-organization process, a number 
of in vitro experiments were performed 0, 0, El IE LJ to 
study interaction of molecular motors and microtubules 
in isolation from other biophysical processes simultane- 
ously occurring in vivo. At large enough concentration 
of molecular motors and microtubules, the latter orga- 
nize in asters and vortices depending on the type and 
concentration of molecular motors. 

In the above experiments the elementary interaction 
processes were identified. After molecular motor binds 
to a microtubule at a random position, it marches along 
it in a fixed direction until it unbinds without appre- 
ciable displacement of microtubules (since the mass of a 
molecular motor is small in comparison with that of the 
microtubule). However, if a molecular motor binds to 
two microtubules (most molecular motors have at least 
two binding sites), it can change their mutual position 
and orientation significantly. In small-scale simulations 
P| , the interaction of rod- like filaments via motor binding 
has been studied, and patterns resembling experimental 
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FIG. 1: a - sketch of motor-mediated two-rod interaction for 
7 = 1/2, b - integration regions Ci,2 for Eq.(|2Jl. 



ones were observed. In Q a phenomenological model for 
the molecular motor density and the microtubule orien- 
tation has been proposed. The model included transport 
of molecular motors along microtubules and alignment 
of microtubules mediated by molecular motors. Simula- 
tions showed that vortices and asters indeed form in this 
model, however, only one large vortex formed in case of 
high density of motors. Ref. generalized this model 
by including separate densities of free and bound molec- 
ular motors, as well as the density of microtubules. This 
model exhibited a transition from asters to vortices as 
the density of molecular motors is increased, in apparent 
disagreement with experimental evidence pj that asters 
give way to vortices with decreasing the molecular mo- 
tors concentration. Somewhat similar approach was em- 
ployed more recently by Sankararaman et al. |10|. again 
with the same apparent disagreement with the experi- 
ment as in Refs. ]^,@. A phenomenological flux- force 
relation for active gels was introduced in . While vor- 
tex and aster solutions were obtained in a certain limit, 
an analysis of that model is difficult because of the large 
number of unknown parameters and fields. 

In Ref. [12| a set of equations for microtubules den- 
sity and orientation was derived from conservation laws 
for microtubules probability distribution function. These 
conservation laws were based on the phenomenological 
expressions for the probability fluxes due to diffusion 
and motor-mediated interactions. The latter, however, 
assumed that tubules are only displaced and rotated in- 
finitesimally in individual interactions, which may not 
be the case in experiments. The model does not produce 
onset of spontaneous orientation for any density of mi- 
crotubules. The authors argue that asters and vortices 
may be created as a result of the "bundling instability" 
|13|. However, this model demonstrated the bundling 
instability at small densities of tubules, and oscillatory 
orientational instability at large densities contrary to ob- 
servations in which the orientational ordering is observed 
at smaller densities and is not oscillatory. In subsequent 
publication [Llj . a derivation of the probability conser- 
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vation equations from microscopic mean-field model of 
forces between tubules and motors was presented, how- 
ever the assumptions made in the course of derivation 
(i.e. infinitely stiff molecular motors) lead to a surprising 
conclusion that filaments do not change their orientation 
during interaction. 

In this paper we present an alternative calculation of 
the microscopic motion of two filaments connected by 
a moving motor in a viscous fluid and show that fila- 
ments do change orientation as a result of the interac- 
tion. In our short publication [l5| we derived a contin- 
uum model for the collective spatio-temporal dynamics of 
microtubules starting with a stochastic microscopic mas- 
ter equation for interacting inelastic polar rods, assuming 
that the density of molecular motors is homogeneous in 
space. Our model differed from the transport equations 
[l2l | in that it treated the interaction between two tubules 
as an instantaneous inelastic "collision" that can change 
the orientation of the filaments significantly. The model 
exhibits an onset of orientational order for large enough 
density of mictorubules and molecular motors jlfij . for- 
mation of vortices and then asters with increase in the 
molecular motors concentration, in a qualitative agree- 
ment with experiment. 

In this paper we present a more detailed derivation 
of our equations starting from a microscopic model of 
tubule-motor interaction. We derive the interaction ker- 
nel from microscopic rules and relate the kernel charac- 
teristics with the properties of the motors. We further 
extend our analysis: we lift the assumption of the ho- 
mogeneous distribution of molecular motors and include 
an additional equation for the evolution of the molecular 
motor density. We also consider the situation when some 
molecular motors are attached to a substrate (usually, a 
glass plate). These modifications allow as to improve an 
agreement with experiment. We explain accumulation of 
molecular motors at the center of an aster by advection 
along microtubules and slow rotation of vortices by the 
interaction of microtubules with the attached motors, as 
it was observed in experiments. 

The structure of the paper is as follows. In Sec. [D] 
we present derivation of the coarse-grained equation for 
orientation and density from the microscopic Maxwell 
model for interacting polar rods. In Sec. Illll we calculate 
the form of the interaction kernel from microscopic rules 
of interaction of two filaments. In Sec. IIVI the effects of 
spatial coupling are considered. First in Sec. IIV Al we 
reduce stochastic equations for the probability function 
to the set of equations for local orientation and density 
of micortubules. In Sec. IIV Bl the stability analysis of the 
isolated vortex and aster solutions is performed and the 
phase diagram of various regimes is presented as the func- 
tion of motor density and anisotropy of the interaction 
kernel. In Sec. [V]we include effects of motors attached 
to the substrate and explain rotation of the vortices and 
onset of large variations of the microtubules density. In 
Sec. I VII we consider effects of variable motor density and 
derive the equation for the evolution of the motors. Tech- 



nical details of derivations are presented in Appendices. 

II. MAXWELL MODEL FOR INELASTIC 
POLAR RODS 

At this stage, molecular motors enter the model im- 
plicitly by specifying the microscopic interaction rules 
between two rods. Since the diffusion of small motors is 
about two order of magnitude higher than that of large 
and heavy mictortubules, in this Section we neglect spa- 
tial variations of the motor density and treat the collision 
rules as spatially homogeneous. Effects of variable mo- 
tor density are considered in Sec. IVII below. Each rod 
is assumed to be of length I and diameter d <^ I, and is 
characterized by the position of its center of mass r and 
orientation angle (f>. 

Consider the orientational dynamics only and ignore 
the locales of interacting rods (an analog of the Maxwell 
model of binary collisions in kinetic theory of gases, see 
e.g. ^3)- We model the motor-mediated inelastic inter- 
action by an instantaneous collision in which two rods 
change their orientations according to the following col- 
lision rule: 

where fa\ 2 are the two rods' orientations before and fa[ 2 
after the collision, and 7 characterizes inelasticity of col- 
lisions. The angle between two rods is reduced after 
the collision by a factor 2j — 1 . Totally elastic collision 
corresponds to 7 = (the rods exchange their angles) 
and a totally inelastic collision corresponds to 7 = 1/2: 
rods acquire identical orientation <f)f 2 = (see 
Fig. nja). Here we assume that two rods only interact 
if the angle between them before collision is less than 
<j>o, \cf>2 — <f>i\ < fa < 7T. Because of 27r-periodicity, we 
have to add the rule of collision between two rods with 
2ir — fa < |02 — 4"i\ < 27r. In this case we have to replace 
fa^ a — > 4>{ a + 7T, fa a — * <p 2 ' a — 7T in Eq. Q. In the fol- 
lowing we will only consider the case of totally inelastic 
rods (7 = 1/2) and fa — tt, the generalization for arbi- 
trary 7 and fa is straightforward, see Appendix 1X1 The 
probability P(4>) obeys the following master equation 

d t P(fa = D r dlP{fa +g f dfadxj> 2 P{fa)P{fa) (2) 

JCi 

x [8{4> - fa/2 - fa/2) - 5(4> ~ fa)] +9 [ dfadfa 

Jc 2 

xP{fa)P{fa)[S{0 - fa/2 - fa/2 - tt) - <S(0 - fa)} 

where g is the "collision rate" proportional to the num- 
ber density of molecular motors m, the diffusion term 
cx D r describes thermal fluctuations of rod orientations, 
and the integration domains Ci, C 2 are shown in Fig^>. 
From the dimensional analysis one finds that the colli- 
sion rate g is of the order of mD r So, since l/D r is the 
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FIG. 2: Stationary solutions P(4>) f° r different p. Inset: the 
stationary value of |r| vs p obtained from the Maxwell model 
© , dashed line - truncated model IKty . 



only time scale in Eq. J2J , and So ~ ^ 2 is the interaction 
cross-section of microtubules. 

Changing variables t — > P r i, P — * gP/D r , w = 4>2~4>i, 
one obtains 



d t P(<« = t^P(0) 



G?U> 



x [P(<j) + w/2)P(<f> - w/2) - P{(j))P{(j) - w)] (3) 

The rescaled number density p = ^ P(4>, t)d(f> now is 
proportional to the density of rods multiplied by the den- 
sity of motors. In the following, an increase of the density 
of molecular motors is reflected in our analysis as an in- 
crease of the number density p. 



A. Orientation instability 

Eq. © possesses uniform steady-state solution P(4>) = 
Pq = p/2tt = const corresponding to isotropic distribu- 
tion of rods. This solution looses its stability with respect 
to anisotropic perturbations with the increase of density 
p. The instability signals the onset of spontaneous ori- 
entation. Substituting solution to Eq. in the form 
P(4>,t) = Pq + £(4>,t), where £ is small perturbation, we 
obtain linear equation for £: 



d t m = dim + f 



m+w/2) 

+ e(0-to/2)-£(0-«j)-e(0))du; (4) 



Looking for the solution to Eq. Q in the form £ ~ 
exp[Afci ± ik(f>], where k ^ is integer, for the growth 
rate A& we find 



(5) 



Afc = p [ — sin(7rfc/2) - 1 



Thus, it follows from Eq. (JSJ that perturbations with 
k = ±1 have the largest growth rate Ai = p(4/7r — 1) — 1. 



The instability ( Ai > ) occurs for the density p > 
p c = 7r/(4 — 7r) w 3.662, and leads to breaking the az- 
imuthal symmetry and formation of anisotropic, i.e. ori- 
ented states. The resulting orientation is determined by 
initial conditions contained in the perturbation £. 



B. Fourier expansion 

Let us consider the Fourier harmonics of the probabil- 
ity density P(<j>): 



-ik<f>\ 



1 



2tt 



— / d0e-***P(^t) 

Z7T 



(6) 



The zeroth harmonic Pq = p/2tt — const, and the real 
and imaginary parts of Pi represent the components t x = 
(cos (j)),T y — (sin</>) of the average orientation vector t, 
t x + ir y — P*. Substituting 10 into Eq.(j2| yields: 



P fc + (k 2 + p)P k = 2ttJ2 Pk- m PmS[TTk/2 - 



(7) 



(here Six) = sinx/x). Due to the rotational diffusion 
term, the magnitudes of high-order harmonics decay ex- 
ponentially with |fc|, see Eq. (JSJ. Neglecting all P& for 
\k\ > 2 we obtain from Eq.® 



Pi 

P 2 + 4P 2 



Pi = P Pi2(4 - tt) - ^P 2 P? 



-P P 2 2tt + 2ttP? 



(8) 
(9) 



Since near the instability threshold the decay rate of 
P2 is much larger than the growth rate of Pi, see Eq. 
(J5J, we can neglect the time derivative P 2 and obtain 
P 2 = API with A = 2ir(p + 4) _1 and arrive at: 



with 



T = €T - A Q \t\ 2 T 



e = piAir- 1 - 1) - 1 » 0.273p - 1 



A = 8A/3 = 



Wit 



3(p + 4) 



(10) 



(11) 



For large enough p > p c = 3.662, an ordering instability 
leads to spontaneous rods alignment. This instability 
saturates at the value determined by p. Close to the 
threshold A « 2.18. 

In order to verify our approximations we solved Eq. 
© for 7 = 1/2, (fio = tt numerically by the finite dif- 
ference method. We find that random initial conditions 
rapidly evolved towards a single-peaked stationary distri- 
bution, the position of the maximum of the distribution 
being determined by initial conditions. Fig. |2 shows 
typical stationary solutions P(4>) obtained from Eq. © 
for different values of p. One sees that P((f>) is weakly 
non-uniform near the critical density and becomes more 
peaked with the increase of the density. From the nu- 
merically obtained distribution we extracted orientation 



amplitude |t| and compared it with the analytical result 
|t| 2 = e/A from Eqs. iJTU jl .lfTT Jl . As seen in the Inset, 
the corresponding values of |r| are consistent with the 
truncated model IjlUfl up to p < 5.5. 

We also studied Eq. © for 4>q < ir. Whereas for <fio 
close to 7r no qualitative difference was found, for smaller 
4>q , e.g. (f>o < 7r/2 we often obtained long- living multi- 
ple peak distributions, with the number of peaks roughly 
7t/0 O - While the distribution possibly relaxes towards a 
single peak, the transient time appears to be very large 
due to exponentially weak interaction between the peaks. 



III. MICROSCOPIC PICTURE OF 
TUBULE-TUBULE INTERACTION 

In order to generalize our model to account for spa- 
tial localization of interaction between tubules, we need 
to introduce a more specific model of interaction be- 
tween two tubules mediated by a molecular motor. 
Namely, we should specify the "collision rules", or a 
relationship between positions and orientations of two 
tubules before and after the interaction via motor attach- 
ment/detachment, and the "collision rate", or the prob- 
ability of the collision to occur given the positions and 
orientations of two tubules. The latter will play a role of 
interaction kernel in the corresponding master equation 
for the tubule probability distribution. 



A. Collision rules 

Here we specify these rules by integrating the equa- 
tions of motion of the two tubules. This calculation is 
based on a number of simplifications. We assume that 
two infinitely rigid rods of equal length I interact with 
one molecular motor. We assume that the motor moves 
with constant speed V along the rods (the results can are 
trivially generalized for the case of V 7^ const) . To sim- 
plify the system even further, we consider a symmetric 
case: the distance of the motor from the center point of 
the rod S, —1/2<S< 1/2 is the same for both rods, see 
Figure [3J Since the size of a motor (ss 30 nm) is much 
smaller than the length of a microtubule (5.. 10 micron), 
we consider a limit of zero motor size. Since the mo- 
tor's bending elasticity is rather small, we approximate 
the motor by a soft spring and prescribe that the force 
F exerted on the tubules due to motor motion is perpen- 
dicular to the bisector of the angle between the tubules 
(i.e. along the motor), which in the symmetric case is 
along the x-axis, see Fig. [3J Even if the symmetry is ini- 
tially broken, and the force is exerted at an angle to the 
x-axis, the force will initiate a relative displacement of 
the tubules in the y direction which will shift the binding 
points in such a way as to restore the symmetry. 

The equations governing evolution of the angle cj> be- 
tween the microtubule and the bisector and the coordi- 
nates X, Y of the center of mass of the microtubule are 
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FIG. 3: Sketch of the interaction of two microtubules 

obtained from balance of torques and forces due to motor 
motion and viscous drag forces 



d t cj> = ^ScaafflF (12) 
d t X = (f jT 4 cos 2 <j> + CI 1 sin 2 4>)F, (13) 

d t Y = (q 1 -^J 1 ) wo. t/> coa </>F (14) 

Here £r,£||,£_L are rotational and translational viscous 
drag coefficients, see Eq. Q55JI below. In the following we 
neglect the anisotropy of the translation friction (£|| = £j_ 
then the equations will simplify considerably 



d t 4> = ^Scos^F 
8 t X = q'F, 

d t Y = 0. 



(15) 
(16) 
(17) 



Additional kinematic equation is obtained from the 
condition that the motor is attached at the distance S 
from the center of tubule, which gives 



X 



-S sin ( 



(18) 



Differentiating Eq. (|18|) with time and using dS/dt = 
V, we exclude F and derive an equation for <f> (note that 
the analysis in this Section is also valid for arbitrary time- 
dependent velocity of the motor V(t) > 0) 



dt 



kVScos((f>) sin(0) 
l + ZcS^cos 2 ^) 



(19) 



where k — £\\/£ r ~ 12/Z 2 , see Eq. lf55f . We make the 
following substitutions: 



r — > kS 2 ; u — > cos(</>) 2 
In new variables Eq (|19f) can be written as 

dr I + tu 
du u(l — u) 



(20) 



(21) 
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Eq. (|21() is linear with respect to r and therefore has an 
exact solution 



C + log(u) 



1 



(22) 



where C is a constant determined by initial conditions. 
Returning to original variables, we obtain 



C + log(cos 2 (</>)) 



kS 2 



(23) 



sin(0) 2 

For small angles <fi Eq. (|23|l simplifies and we obtain 



VI + kS 2 



(24) 



where 0o an d <So are the initial conditions at t = 0. 

For the final angle obtained when the motor reaches 
the end of the microtubules (S = 1/2) we obtain 



^1 + kP/4 



(25) 



As one sees from Eq. I|25|) . the final angle (/> depends 
on the initial angle </>o and the initial attachment posi- 
tion So- Assuming that the probability of attachment of 
the motor is independent of the position along the micro- 
tubule S, in the small angle approximation for average 
angle 



(4>) = I' 



1/2 



4>(S )dS 



we obtain 



(4>) 



1/2 



asinh(Vfcl2/2) 



(26) 



VWyJl + kP/A 
Thus, the averaged change in the angle e = (</>)/ '0 O is 



(27) 



_ 1 asinh(v / M2/2) 

e ~ 2 Vw^/TTWjt 



(28) 



Obviously, for fcZ 2 — > oo the relative angle change e — > 
1/2. Correspondingly, the restitution coefficient 7 = 
(1 — e)/2 — * 1/4, in agreement with our assumptions on 
inelastic collision between the rods. For the case of rigid 
rods we obtain from Eq. (|55|l that kl 2 w 12, which gives 
e « 0.67, or 7 ss 0.17, which is considerably smaller than 
fully- inelastic case of 7 = 1/2. With the increase of k the 
coefficient 7 increase, e.g. for k = 30 one obtains 7 = 0.2. 

For arbitrary <j>$, in order to find the average angle 
change we solved Eq. , J5BJ| numerically (see Fig. 0J. 
There is a very week dependence of 7 on the initial angle 
4>. In a wide range of angles Aip < 0.75-7T, parameter 
7 « 0.2 and then 7^0 for — > tt. 
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FIG. 4: The effective restitution coefficient 7 calculated from 
Eq. (12(il for k = 12 (solid line). For comparison 7 at A: = 30 is 
also shown (dashed line). Dot-dashed line indicates the limit 
of fully inelastic interaction 7 = 0.5. 



B. Collision rate 

Now we turn to the calculation of the collision rate 
IF(ri, r2, 0i, <t>2) between two tubules with center-of- 
mass positions iq,r2 and orientation angles 4>\,4>2- We 
consider a translationally and rotationally invariant sys- 
tem, so the collision rate depends only on the position 
and orientation differences, W(r, <f>), where r = 17 r2, 
and, as before, <j> = (4>i — </>2)/2. 

Microtubules interact via a motor attachment only if 
they intersect, which is expressed by the condition 



ki - r 2 | < ^sin< 



This gives the collision rate in the form 

I^ = 5 e(|r 1 -r 2 |-Z/2sin0) 



(29) 



(30) 



where is the 0-function, and g is the probability of a 
motor to bind to both tubules given that they intersect. 
The latter is directly proportional to the concentration 
of molecular motors bound to a tubule. If the concentra- 
tion of the motors along the tubules were uniform, the 
collision rate would be uniform inside the parallelogram 
defined by Ea. (|30|l . However, due to the transport of 
motors along the tubules, their concentrations increase 
towards the polar ends of the filaments. 

In order to compute the inhomogeneous motor distri- 
bution along a filament, we assume that the motors can 
be in two different states, bound and free. The concen- 
tration of free motors nif is a function of the coordinate 
along the tubule S, (—1/2 < S < 1/2), perpendicular 
coordinate r±_ (we consider a two-dimensional domain), 
and time t. Bound motors are localized on the tubule 
itself, so their concentration is Mf,(S, t)S(r±). 
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The inhomogeneous motor distribution can be evalu- 
ated from the following microscopic equations for motor 
attachment/detachment and advection processes: 

d t irif = (p ffM b ~p on m f )d(r^) + DV 2 mf (31) 
d t M b = -p off M b +p on m f (S,0,t)-d s {VM b ) 

These equations, formulated in terms of the concen- 
tration of bound/free motors M b ,mf, describe random 
binding/unbinding of the motors with the probabilities 
Pon,off] diffusion of free motors (diffusion coefficient D), 
and the advection of bound motors with the velocity V 
along the tubule. The parameter p ff characterizes pro- 
cessivity of molecular motors: large p ff corresponds to 
small processivity, motor unbinds soon after it binds to 
a filament. 

According to experiments (see e.g. 0)j while multi- 
ple motors can be attached to a single tubule, only one 
motor can be attached per elementary binding site (which 
represents a section of approximately Iq = 10 nm along 
the tubule) . This leads to a kind of hard-core repulsion 
which in the simplest approximation can be taken into 
account by introducing local "pressure" P of bound mo- 
tors, and modifying the transport speed V , 



v = v - v d s p 



(32) 



where r/ is an effective mobility. Pressure P diverges as as 
the bound motor density approaches densely packed limit 
of one motor per binding site M b — > Mo ~ 1/Iq- We will 
adopt a simple generic expression for the pressure as a 
function of the bound motor density (see for comparison 
expressions for pressure in granular hydrodynamics near 
closed-packed density l2ll E^l 



P 



M h T 



1 - M b /Mo 



(33) 



The "temperature" T is determined by fluctuations of 
bound motors on the tubule and is typically small, so 
the pressure can be neglected everywhere except where 
the density is very close to the dense limit. 

In the stationary state Eqs. I|32[l assume the form 

(p ffM b - p on m f )S( ri _) + DV 2 m, f = (34) 
-PoffM b + p on m f {S, 0) = d s [V M b - V M b d s P{35) 

Since the diffusion constant D of free motors is large, 
we can neglect the inhomogeneity in the free motor dis- 
tribution and assume m/ = const in Ea . 1)35(1 . Eq. 1(35(1 
has to be solved with the boundary condition M b = 
at S = -1/2. At the end of the tubule S = 1/2 the 
"exit" flux of bound motors VM b is determined by the 
detachment probability p en d of the motor, resulting in 
the condition VM b = p end loM b . For small T, the den- 
sity of bound motors has two distinct phases: low-density 
"gas phase" near the beginning of the tubule (S = —1/2), 
and a high density "solid phase" near the end (S — 1/2) 
[23| . The location of the boundary between these phases 



can be found by equating the fluxes of bound motors in 
the two phases. In the low density phase, the pressure 
term can be neglected due to small T, and the solution 
has the form 



Mb = 



P on m s 

Poff 



1 — exp 



Poff 

V 



(S + l/2) 



(36) 



Typically, PoffV^l 3> 1, so the density saturates very 
quickly to the equilibrium value M e — p ori p~jjnif. This 
solution corresponds to a constant flux of motors along 
the tubule, 



F l = V M e = 



Poff 



(37) 



In the solid phase at a very low temperature T, the 
motor density is very close to Mq. Thus, at the end of 
the tubule, the flux of motors is equal to the number of 
motors leaving the tubule in a unit time p end Mo- Ac- 
cording to Ea. (|35|l . the flux of motors in the solid phase 
is a linear function of the coordinate, 

F 2 = PendhMo + [p f fM - p on m f ] {1/2 - S) (38) 

The two phases are separated by a narrow interface (the 
width of the interface vanishes when T — > 0) whose po- 
sition iSo is determined by equating these two fluxes, 
Fi = F 2 , 



VoPonmf 
Poff 



■ p end l M + \p ffM -p on m f }(l/2-So) (39) 



This yields the following expression for the interface po- 
sition 



So 



1 p end l M - V p on p of 1 f mf 

2 PoffM - PonTUf 



(40) 



Obviously, S grows with p end , and at p end = 
VoPonmflpoffloMo} 1 , we obtain So — 1/2, i.e. the solid 
phase disappears. 

Thus, the bound motor density M b is approximately 
described by the following step function 



M b (S) w M e + (M - M e )Q(S - So) 



(41) 



The inhomogeneous distribution of bound motors di- 
rectly leads to the anisotropy of the collision rate. The 
coordinates along the microtubules Si, 2 are related to the 
positions ri ^ of the center of microtubules as follows 



Si, 2 = ni )2 • (r - ri, 2 ) 



(42) 



The collision rate is proportional to the sum of the bound 
motor concentration times the cross-linking rate go, 



g = <7 [A/ 6 (S?) + M b (S° 2 )} 



(43) 



7 



where S® 2 are the values of S% t 2 at the intersection point. 
Excluding r from Eq. (|42|) one obtains 



S? 
S 2 ° 



(£1 - r 2 )(n 2 - ni(nin 2 )) 
1 - (nin 2 ) 2 
(ri ~ rgKni - n 2 (n 1 n 2 )) 
1 - (nin 2 ) 2 



(44) 



Non-uniform Mb(S) produce anisotropy: 
5(ri,ni,r 2 ,n 2 ) ^ g(n, n 2 , r 2 , ni). However, the 
collision rate Eq. i|43|) with the step-wise expressions 



for S® 2 <|41[) is awkward and impractical for further 
calculations. In the subsequent section we will not use 
the exact form (|3UI) with Ea. l|43|) as a kernel in the 
master equation, but replace it with a more simple form 
which nevertheless retains the main features of H30|) . 143f) . 
namely, localization and anisotropy, 

W « Wo (n - r 2 ) (1 - /^(ri - r 2 )(m - n 2 )) (45) 

where the symmetric part of the kernel Wq is of the Gaus- 
sian form 



VF (r) = ^cxp[-r 2 /& 2 



(46) 



with the spatial scale b w 1/2. The dimensionless param- 
eter (3 characterized the collision rate anisotropy. The 
interaction kernel in this from was proposed by us on the 
symmetry grounds in Ref. |T^ . 

While the form (|45|l cannot be rigorously derived from 
(|30|1 . (|43|) . the anisotropy coefficient (3 as a function of 
kinetic parameters can be estimated from the expression 
(I41|l . First of all, we approximate the step function in 
Eq. (|41() by the linear function 



M h {S) 7zM + aS 



(47) 



where the mean density M and mean slope a are calcu- 
lated by the least mean square method using Eq. (|41|) . 



" 2 MbdS = m °+ m - 

4/2 2 
1/2 



(M -M e )^ 



p 



1/2 



p 



--5 2 



(48) 



To evaluate the effective kernel we substitute Eq. Q47JI 
into Eq. Ij43(l and using relations (|44|l . obtain 

g » .go [2M + a(5° + S 2 )] 

= go (2M-a ^-^-^ ) (49) 
V 1 - nin 2 y 

As one sees from Eq. (|49fl , it coincides with phenomeno- 
logical kernel Eq. I|45|l up to the factor (1 — nin 2 ) 
in denominator. The value of the dimensionless kernel 
anisotropy is then 







ctl 
2M 



(50) 



Assuming that the density in the solid phase Mq is 
much larger than the density in the gas phase M e , 
fi = Mo/Me ^> 1, we obtain the following estimate for 
the anisotropy (3 for p en d —* 0: 



/3« 3 (--^) =3 Fo 7 Pe "^° . (51) 
2 I J IPoffH 



The anisotropy is maximum for p e nd = 0, {(3 « 
Vo/pofflfx), decreases with the increase of p en d and van- 



ishes (at this approximation) at p e 



Vo/filo- 



We can estimate the parameter (3 for different type of 
motors using the data from Ref. 0, . The parameters 
for kinesin and NCD are: V = 1 /i/sec, p on = 20 sec -1 , 
Po// = 0.5 sec -1 , and p en d = 70 sec -1 for kinesin and 
Pend — 2.5 sec -1 for NCd. Projectedftwo-dimensional) 
density of free motors to/ in Refs. 0,H| was taken m/ = 
0.05 — 2 /im -2 . For the linear density of bound motors 
we obtain M e = Pon'mfd^/poff and Mq = 1/Zo, where 
do ~ 2/o is the diameter of microtubule. For parameter /i 
obtain [i = M /M e = p Q f f I {d l p on mf) ~ 10. .500 ^> 1. 
Thus, for NCD-like motors when lo[ip en d < 1 we obtain 
the anisotropy parameter (3 (depending on the density 
ratios /i) 



(3ncd 



3V 



10 -;i ..10" 



(52) 



It follows from Eq. I|52|) that the kernel anisotropy in- 
creases with the increase of the concentration of free 
motors. Correspondingly, since for kinesin parameters 
Pend > Vop on mf[p f /Zo^o] 1 an d no solid phase is 
formed, the anisotropy coefficient (3 is essentially zero. 

In this Section we considered only one mechanism con- 
tributing to the kernel anisotropy: inhomogeneous distri- 
bution of bound motors. Possibly there are other mech- 
anisms affecting the anisotropy of the interaction, for ex- 
ample, finite bending rigidity of the microtubules may 
also contribute to both the collision rules and the col- 
lision rate. However we leave this interesting topic to 
further studies. 



IV. SPATIAL LOCALIZATION OF 
TUBULE-TUBULE INTERACTION 



To describe the spatial localization of tubule-tubule 
interaction we introduce the probability distribution 
P(r, (f>, t) to find a rod with orientation </> at location r at 
time t. The master equation for P(r, <fi, t) can be written 
as 
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d t P{r^) =dlP(r,<f>)+d i D ij d j P(r,<t>)+ / / dndr 2 / dw [W^,^^ + w/2, (j>- w/2) 



xP(r 1 ,4> + w/2)P(r 2 ,<t>-w/2)5 ( ^ t - - r 



-Wfa.ra, ^, - w)P(r 2 , 0)P(r 1; - w)<5 (r 2 - r)] (53) 



r 



where we performed the same rescaling as in Eq. © using 
go instead of g and dropped argument t for brevity. We 
normalized the probability P — ► P/l 2 . 

Unlike the Maxwell model equation J3J), Eq. I|53[l con- 
tains two diffusion terms (translational and angular), and 
the motor-mediated tubule-tubule collision integral now 
contains an interaction kernel which depends on the rela- 
tive tubule positions and orientations. The angular diffu- 
sion coefficient D r and the translational diffusion tensor 
Dij are known from the polymer physics 24] : 



Dij = — (D\\mnj + D±(6i 



D 



D± = 

D r = 



k B T 
k B T 
4fcgT 



(54) 



where £ii,£_L,£r are corresponding viscous drag coeffi- 
cients. For rigid rod-like molecules, 



2irr] s l 

foiWd) 



a =2^||; & 



nr] s l 3 



31og(//d) 



(55) 



where r] s is shear viscosity . 
the rotational diffusion time t 



Since we scaled time by 
-> D r t, in new variables 



the translational diffusions assume a very simple form: 
D|| = 1/24, Di. = 1/48, see Eq. Note that the 

drag coefficients are slightly modified for thin films and 
membranes, see e.g plf . 

The last term of Ea. l|53|l describes motors- mediated 
interaction of rods. While the results of the previous sec- 
tion indicate that the angular "inelasticity coefficient" 
7 is less than 1/2, we postulate here that after the in- 
teraction, the two rods acquire the same orientation, 
4> = {4>i + ^2 ) and the same spatial location in the mid- 
dle of their original locations, r = (ri + r2)/2. These 
assumptions are made to simplify the calculations and 
final equations, however generalization to arbitrary colli- 
sion inelasticity is straightforward (see Appendix |A"|) . 



A. Continuum equations 



After integration over the S functions Eq. (|53() assumes 
the form 



dP d 2 P 



dt d(j) 2 
where nonlinear terms 



diDijdjP + Z + pZi (56) 



and 



Z$= dn / dti;[2Wb(2(n-r))P(r 1 ,0 + t«7 s t)P(2r-n,0 + i£;(7-l),«) 

J J —<j>0 

-W {v 1 -v)P{v,4>,t)P{r u 4>-w 1 t)] 



(57) 



Zi= I dr x I ' ° dw [2W (2( ri - r))(n - r) • (m - n)P(n, <j> + vyy, i)P(2r -r u <j> + w( 7 - 1), t) 

J J —<t>o 

-W (ri - r)(n - r) • (n x - n)P(r, </>, t)P(r u <j> - w, t)] (58) 



are generated by the collision integral in Eq. (|53|l . 

To proceed, we again perform the Fourier expansion of 
the probability distributions in ip and truncate the series 
at \n\ > 2. Now 27rPo gives the local number density 



p(r,t), and P±\ the local orientation r(r, t). The inte- 
gration of the diffusion term in Ea. (|56|l generates linear 
terms, and the nonlinear terms Z$, Z\ (see Appendix IE|) 
produce nonlinear terms in the corresponding equations 
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for p, t. 

After rescaling space by I, and introducing dimension- 
less parameters B = b/l characterizing the width of the 
interaction kernel and H — j3B 2 characterizing normal- 



ized strength of anisotropy of interaction, we arrive at 
the set of equations for coarse-grained local density p 
and orientation r 



dtp 



£_ 

32 



B 2 p 2 
16 



irB 2 H 



1G 



[3V • (rV 2 p - P V 2 t) + 2d t (djpdjn - d t pd 3 T 3 )} 



7 P pB 4 
256 



t r = ^r + iv(V.r) 



er - A \t\ 2 t - H 



Vp 2 



16tt 



- t(V-t)--(tV)t 



An 



(59) 
(60) 



These equations generalize Ea. (|10(l for the case of spa- 
tially localized coupling. For simplicity the last two terms 
in Eqs. (|59|I . (|6U|) have been linearized near the mean den- 
sity po = (p), otherwise more complicated expressions 
given in Appendix [B] are needed. This approximation 
is justified by the fact that in the relevant range of pa- 
rameters B, H the density variations are small compared 
to the variations of the orientation r. The last term 
in Eq. H59|) regularizes the short-wave instability when 
the diffusion term changes sign for pa > Pb = l/AB 2 . 
This instability leads to strong density variations asso- 
ciated with formation of dense microtubule bundles (see 
Figs. El 03 below) which is also observed experimentally 
for large density of molecular motors. 



B. Stability of Asters and Vortices 



If B 2 H -c 1, the density modulations are rather small, 
and Eq. (|60|l for orientation r decouples from Eq. (|59|l . 
It is convenient to rewrite Eq. (|60|l for complex variable 

■0 = T x + iT y : 

dti/> = eip- A |^|V + £>iV 2 i/> + £> 2 V"V 



+H I 1 7T - -J VReVV>* + -(Re^*V)^J (61) 

where operator V = d x +id y , D\ = 1/32+poB 2 /An, D 2 = 
1/192. Eq. (|61|l is similar to generalized Ginzburg- 
Landau equation known in the context of superconduc- 
tivity, superfluidity, nonlinear optics, and pattern forma- 
tion, see e.g. [27j • Let us focus on the structure and the 
dynamics of radially-symmetric solutions of l|61|l which 
can be sought in polar coordinates r, 9 in the following 
generic form: 



ip = y/A /eA(r) exp(i6) 



(62) 



where the complex amplitude A(r) = $(r) exp[iip(r)], 
and the phase <p(r) is a real function. The solution 
tp(r) = 0, 7r corresponds to asters, and the solution with 
ip(r) 7^ describes vortices, see Fig. |SJ Transitions be- 
tween asters and vortices can be examined in the frame- 



work of a one-dimensional problem for the complex vari- 
able A(r), 

d t A = DiA r A + D 2 A r A* + (l - \A\ 2 ) A 
+H ( ai^ReV r A + a 2 d r AReA + — ) (63) 



with the following differential operators 



A,. = d 2 + r^dr - r~ 2 ; V r = d r + r' 1 (64) 



and parameters 



Oi 



(tt-8/3)/, 



0.321, a 2 = 8/3a 



1.81 



Here we rescaled time t — ► t/e and space variable by 
r — > rj^fl. 

The aster and vortex solutions obtained by numerical 
integration of Eq. (|63|l for certain parameter values are 
shown in Fig. [7| Vortices exist only for small values of 
H and give way to asters for larger H or larger p. For 
H = 0, Eq.jnSI) reduces to a form which was studied in 
ptjf . It was shown in j2|| that the term A r A* favors 
vortex solution (ip — ±7r/2). In contrast, terms propor- 
tional to H select asters with ip = tt (aster for ip = is 
unstable for the anisotropy parameter H > 0). Note that 
ip = 7r corresponds to asters with the direction of arrows 
towards the center, as it is shown in Fig. [S] Since we as- 
sociated the direction of the vector r with the direction 
of motion of molecular motors along the microtubules, 
the aster with ip — w corresponds to the experimental 
situation: motor moves towards the center of the aster. 
Increasing H leads to gradual reduction of <p, and at a 
finite Hq(pq), 4>{r) = tt, i.e. the transition from vortices 
to asters occurs. For < H < H , the vortex solution 
has a non-trivial structure. As seen in Fig. [7| the phase 
<p — > 7r for r — > 00, i.e. vortices and asters become indis- 
tinguishable far away from the core. 

The transitions between asters and vortices can be 
studied in the framework of linearized Eq. I|63|) . For 
this purpose the solution to Eq. (|63|l can be sought in 
the form 



A(r) = <f>(r) + iw(r) exp(At) 



(65) 
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FIG. 5: Schematic representation of orientation fields r for 
three different values of ip; aster [ip — 7r); generic vortex 
(•7t/2 < ^ < 7r) and ideal vortex (tp = ±7r/2). 
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where small real perturbation w obeys a linear equation 
L = Xw with operator 

L = DA r + (1 - $ 2 + aiiJV r $) + a 2 iJ$V r (66) 

(D = Di — D 2 ) with zero boundary conditions at r = 
0, 00. This eigenvalue problem can be solved by the 
matching-shooting method. A positive eigenvalue A cor- 
responds to the emergence of a non-zero phase (j>(r), i.e. 
a vortex. 

The resulting phase diagram of the continuum model 
(I59|l . (|60|l is shown in Fig. HO The solid line i?o(po) sepa- 
rating vortices from asters is obtained from the solution 
of the linearized Eq. Ij63(l by tracking the most unstable 
eigenvalue A of the aster. The dashed line corresponds 
to the onset of the orientation instability, po = Pc- The 
lines meet at the critical point H c — Ho(p c ) above which 
vortices are unstable for arbitrary small e > 0. The phase 
diagram is qualitatively consistent with experiments, see 
Ref. [Bj|: for low value of kernel anisotropy H < H c 
(which according to our estimates in Sec. IIII Bl corre- 
spond to kincsin-like motors with very small anisotropy 
value (3) increase of the density po first leads to formation 
of vortices and then asters. For H > H c (which appar- 
ently correspond to Ned motors with large anisotropy) 
only asters are observed. For large density values in addi- 
tion to orientation instability one observes density insta- 
bility po > pb — 1/4B 2 produced by Eq. J^SJl when the 
diffusion term changes sign. Numerical studies of the full 
system (15911 , l|60f) indeed indicate formation of extended 
domains of high density not associated with the asters, 
see Fig. HO While our model ((23), |S3> yields the den- 
sity instability and bundle formation in accordance with 
experiment, we anticipate only qualitative agreement in 
this regime because the model itself is derived in the low 
density limit when only binary interactions are included. 

C. Interaction of asters and vortices 

For H 7^ well-separated vortices and asters exhibit 
exponentially weak interaction. For asters it follows from 
the fact that L is not a self-adjoint operator. To investi- 
gate the interaction between asters we need to examine 
the asymptotic null-space of the adjoint operator Iy for 
r — > 00 (see for details 27]). After simple algebra we 



FIG. 6: Phase boundaries obtained form the linear stability 
analysis of the aster solution for B 2 — 0.05, dashed line shows 
bundling instability limit po — pt = 5. Inset: Position of 
cr 
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FIG. 7: Stationary vortex and aster solutions r x + ir y — 
$(r) exp[i6> + i<p(rj] to Eq. for p = 4, B 2 = 0.05. 



obtain that for r — > 00 the adjoint operator is of the 
form 

D = Ddl - a 2 Hd r (67) 

Substituting the solution to Eq. l|67|l in the form ~ 
exp[pr] we obtain that there are two solutions: p = 
which describes the neutral translation mode, and the 
non-trivial solution 

p = a 2 H/D (68) 

Since due to interaction the well-separated asters produce 
only small perturbations to their shape, these perturba- 
tions can be treated in linear approximation, and the 
exponent l|68l) characterizes asymptotic screening of the 
interaction between the asters analogous to the interac- 
tion of sp iral waves in the Ginzburg-Landau equation, see 
Ref. [23 for details of analysis. Thus we obtain that per- 
turbations produced due to interaction of well-separated 
asters decay as w ~ exp[— r/Lo] with the screening 
length L = l/p, or in original units Lq = D/a 2 H^/e. = 
D/a 2 H^J p/p c — 1 (see for details 27J). Screening length 
Lq diverges for H — > and at the threshold po — > p c . 
Similar analysis can be performed for vortices. 
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FIG. 8: Composite image of the density (colors) and orienta- 
tion (arrows) fields in the regime of density instability. Den- 
sity changes from p ma x ~ 10 (red) to p m in ~ 4 (dark blue). 
Parameters: B 2 — 0.05, po — 6, H = 0.125, domain of inte- 
gration 80 x 80 units. 




FIG. 9: Orientation r for vortices (H — 0.006, left) and asters 
(H = 0.125, right) obtained from Eqs. 1591 . 16UI . Color code 
indicates the intensity of |t| (red corresponds to maximum 
and blue to zero), B 2 = 0.05, po = 4, domain of integration 
80 x 80 units, time of integration 1000 units. 



D. Numerical solution of full system 

We also studied the full system (|59~|) . i jrJDIl numerically. 
Integration was performed in a two-dimensional square 
domain with periodic boundary conditions by the quasi- 
spectral method. For small H we observed vortices and 
for larger H asters, in agreement with the above anal- 
ysis. Since vortices exist for smaller values of H, their 
screening length L is larger than for the asters because 
L ~ l/H, see Eq. (|68|) . Thus, the vortices interact 
stronger and are more keen to annihilate than asters. As 
seen in Fig. |5| asters have a unique orientation of the 
microtubules (here, towards the center). Asters with the 
opposite orientation of r are unstable. 



In large domains asters form a disordered network of 
cells with the cell size of the order of Lq. Neighboring 
cells are separated by the "shock lines" terminated by 
saddle-type defects. The pattern of asters resembles a 
"frozen" glass state of spirals observed in the complex 
Ginburg-Landau model [23, [28j . Starting from a random 
initial condition we observed initial merging and annihi- 
lation of asters. Eventually, annihilation slows down due 
to exponential weakening of the interaction of asters. For 
the same integration time the number of vortices is typ- 
ically smaller than the number of aster due to the fact 
that the screening length of asters is smaller. 



E. Drift instability of the asters 

In experiments 0, || asters often are not stationary: 
they drift and coalesce. Surprisingly, in our numerical 
investigations of Eqs. I|59|l . I|60|l we also observed that 
typically the center of an aster is unstable and develops 
a spontaneous acceleration instability, see Fig. ^| This 
instability is reminiscent of the instability of the core of 
spiral waves in the complex Ginzburg-Landau equation 
in a large dispersion limit, see Ref. [23. This instability 
associated in Ref. [2^ with the exponential growth of lo- 
calized mode in the form w\(r) exp[i#], which is similar to 
the translation mode and results in the displacement of 
the core. We have found that the instability can be sup- 
pressed by increasing the coefficients in front of the last 
term (~ V 4 p) in Eq. H59fl . In this context it is possible 
that the acceleration instability is just an artifact of the 
approximations made in the course of derivation, such as 
binary character of interactions of microtubules, small- 
gradient expansion, etc. Furthermore, for high density 
of the microtubules the prefactor in form of the cut-off 
term ~ V 4 p to should grow rapidly and dampen the in- 
stability. While at the moment there is no unambigous 
experimental evidence of this instability (e.g. the drift 
of asters could be also attributed to gradients of micro- 
tubules or motor distributions, effects of the boundaries 
etc), we believe that this instability can be found in cer- 
tain experimental conditions. To study the drift insta- 
bility we prepared by the initial condition axisymmetric 
aster solution perturbed by the small amplitude noise. In 
the course of motion the solution breaks the axial sym- 
metry, typical structure of the moving aster is shown in 
Figure EH There is a small but noticeable (about 10%) 
increase of the density p and the amplitude of orientation 
\t\ behind the aster, for the immobile solution the posi- 
tion of the zero of r and maximum of p coincides. The 
instability accelerates collisions and coalescence of asters. 
However the growth rate A of the instability appears to 
be very small and aster solutions are well-preserved for a 
very long time (several hundreds of dimensionless units). 
Fig. ^Jshows the velocity of aster core vs time for the pa- 
rameters of Fig. I1UI One clearly sees initial exponential 
growth of the aster velocity. 
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FIG. 10: Drift instability of asters for H = 0.125, B = 0.06, 
po = 4, size of the image 40 x 40 at the moment of time 
t = 300. Left image shows \t\, right image shows p, arrow 
indicates the direction of drift. Color code: blue corresponds 
to zero, red corresponds to maximum. 




FIG. 11: Aster core velocity vs time for parameters of Fig. 
1101 Solid line shows v y , dotted line shows v x , dashed lines 
depict exponential fit v ~ exp(Xt) for first 150 units of time. 



V. EFFECTS OF MOTORS ATTACHED TO 
THE BOTTOM PLATE 

In the previous Sections we considered microtubules 
interacting with molecular motors freely floating in the 
solvent. However, in in vitro experiments it is difficult 
to prevent attachment of some fraction of motors to the 
bottom of the cell with one of their two heads. The other 
(free) head of the attached (absorbed) molecular motor 
then may bind to a microtubule and push it in the di- 
rection opposite its orientation. This effect was observed 
experimentally in Ref. |30l | (refereed to as microtubule 
gliding assays). 

The effect of attached motors can be easily incorpo- 
rated in the master equation 



dP 
~&t 



d 2 P 



0,1 ),,(), I' 



aV(nP) + Z a + Zi(69) 



in the direction opposite to their orientation vector n = 
[cos </>, sin </>] . Terms Zo,i remain unchanged. It is easy to 
check that the drift term will generate additional linear 
terms airVr in Eq. and a(47r) _1 Vp in Ea. (|6Tj|) . 

It is useful to write Eq. (|59|l in the form of the mass 
conservation law 



Pi 



-VJ 



(70) 



Here a is the fraction of the attached motors, and the 
term nVP account for the transport of the microtubules 



with the corresponding mass flux J. The anisotropic 
part of the kernel generates mass flux which is second 
order in the gradients of p, r. The the lowest-order term 
J ~ r in the expression for flux is generated by the mo- 
tors attached to the substrate and is similar to that of 
self-propelled particles, see e.g. J2(|. In the situation 
considered in Sec. IIV Al this term is prohibited by the 
momentum conservation: the molecular motors produce 
only internal forces which cannot displace the center of 
mass of the system. However this is not the case when 
some of the motors are attached to the substrate. In Ref. 
[hij similar contribution to the flux were attributed to the 
net displacement of the center of the microtubule pair due 
to the anisotropy of the viscous drag coefficient. However 
this pure hydrodynamic effect is probably smaller than 
the advection produced by the motors absorbed at the 
substrate. 

While the fraction of absorbed motors a might be 
small, it still can produce a considerable effect because 
it generates the lowest-order transport term in Eqs. 
l(59 )l .l(50 )) . Numerical studies of Eqs. l(59 )l .l(50 )l with ad- 
ditional a-terms reveal that the qualitative features are 
not very sensitive to the presence of these terms for small 
a <C 1, as long as the diffusive transport in the equa- 
tion for the density p dominates advection. However, for 
moderate a we observed that the aster and vortices be- 
come even less localized, see Fig. This derealization 
is due to the fact that the absorbed motors advect the 
mocrotubules in the direction opposite their orientation. 
Consequently, these motors move the tubules from the 
asters and make a small depression of density for a =/= 
contrary to the density peak for a = 0, compare images 
on Fig. Similar results are also obtained for vortices. 
Remarkably, the suppression of density of microtubules 
in the core of vortex is observed experimentally, see Fig. 
2a in Ref. @. 

The absorbed motors, resulting in the displacement 
of the center of mass of the microtubules system, may 
explain rotation of vortices absent in our previous anal- 
ysis. Indeed, since these motors generate net motion of 
individual microtubules with the velocity ~ a, they can 
support rotating configurations similar to that observed 
in the system of vibrated rods j2(| . Obviously no rota- 
tion anticipated for asters due to pure radial orientation 
of microtubules: the forces induced by motors attached 
to substrate will be compensated by "pressure" gradi- 
ent due to redistribution of density of microtubules. In 
contrary, the rotation is present for vortices. Far away 
from the core the distinction between vortex and aster 
disappears, the rotation is localized only at the core of 
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FIG. 12: Comparison between density distributions (colors) 
without (upper image) and with drift term (a — 0.004). Ar- 
rows show corresponding orientation, overall change in the 
density about 5%, B 2 = 0.05, po = 4, H = 0.1, domain of 
integration 80 x 80. 



the vortex where the phase <fi is different from tt. Since 
the amplitude of orientation vector r grows almost lin- 
early from the vortex core and reaches asymptotic value 
To ~ yJe/Aq at the distance about 1 — 2 dimensionless 
units (see Fig. Q|, the rotation frequency of the vortex 
core u olt q . Indeed, rotation of the vortex core was 
observed experimentally. For the parameters of our nu- 
merical studies the frequency lu is very small due to the 
smallness of a = 0.004, thus during the time of numeri- 
cal experiment (< 1000 dimensionless units of time) the 
vortex core turned only the fraction of full circle. 



VI. INHOMOGENEOUS DISTRIBUTION OF 
MOTORS 

In previous Sections we always assumed a homoge- 
neous bulk distribution of molecular motors (however we 
took into account local inhomogeneity of bound motor 
concentration on the scale of a single tubule to account 
for the collision rate anisotropy in Sec lIII Bl) . This as- 
sumption was justified by the fact that the diffusion of 
motors is about two orders of magnitude larger than of 
microtubules. However experiments indicate that even 
despite this strong diffusion, molecular motors aggregate 
in the core regions of asters and vortices due to the di- 
rected transport of motors by microtubules . 

To describe the dynamics of the motor concentration 
we again invoke the equations for free m/ and bound 
m b motor concentrations, but unlike Sec. IIII Bl wc 
will coarse-grain these distributions on the scale much 
larger than the size of individual filament similar to Ref. 
(see also H, 0, 0) The populations m&, m/ obey 
the advection-diffusion equations jjj (compare with Sec. 
ITTTBTi 

d t m f = DV 2 m f ~ p(p on m f -p off m b ) 

d t m b = -(Vm b T + p(p on m f -p off m b ) (71) 

where p on ,p°ff are the rates of binding/unbinding of mo- 
tor to the microtubules, D, £ are diffusion/advection co- 
efficients accordingly. 

If we assume that the distributions of rrif , m b are 
smooth and the binding/unbinding rates are large, then 
the r.h.s. Eqs. I|71|l is dominated by the last term de- 
scribing binding/unbinding of the motors, leading to the 
local balance relation between m f , m b 

p on m f « p off m b (72) 

Then we can reduce system (|71|l to a single equation for 
the total motor density m = rrif + m b : 

d t m = D V 2 m - CoVtot (73) 

where D = Dp°^/(p°^ + p on ) and Co = (p on /(p off + 
p on ). 

Accordingly, we need to modify the expression for the 
interaction kernel Eq. I|45|l in order to include the effect of 
the motors. The simplest way to include inhomogeneous 
motor density into the kernel is the following, 

W m = m W{ ri - T2 ,<h. - (74) 

where W is given by Eq. (|45|l . Taking the motor con- 
centration in the middle point (ri + r2)/2 is necessary 
to preserve the mass conservation law. Repeating the 
calculations presented in the previous sections one can 
derive equations similar to Eqs. H59|M60fl but with the 
motor density as an independent field. However the re- 
sulting equations are very cumbersome, especially for the 
transport term in Eq. I|59|) . 
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One can simplify the problem considerably utilizing 
again the fact that the motor diffusion is high, and there- 
fore, the distribution of m is smooth. Then, one can ne- 



d t p = 



d t T = 



Thus, motor density is included in the lowest order 
in gradient expansion. We also included terms ~ a de- 
scribing the transport of microtubules by the absorbed 
motors. Again for simplicity we replaced the motor den- 
sity m by its mean value in the last terms in Eqs. 

We carried out numerical studies of Eqs. J75J) , JJJJJ . 

The values of the parameters -Doi£o can be estimated 
from the experimental conditions, in our dimensionlcss 
units Do ~ 1..5 and £o ~ !• 

Selected results are shown in Fig. 1131141 In agree- 
ment with experiment, we observed that motors tend 
to accumulate in the center of an aster or a vortex, see 
Fig ll3ll4l Otherwise, qualitative behavior of formation 
of asters and vortices remains the same. As seen in Fig- 
ure El initial multi- aster state coarsens and leads to the 
formation of a network of large asters separated by the 
domain walls. 



VII. CONCLUSIONS 

In this paper we derived continuous equations for the 
evolution of microtubule concentration and orientation 
due to their interaction via molecular motors. We found 
that an initially disordered system exhibits an order- 
ing instability qualitatively similar to the nematic phase 
transition in ordinary polymers at high density. The im- 
portant difference is that here the ordering instability is 
mediated by molecular motors and can occur at arbitrary 
low densities of microtubules. At the nonlinear stage, 
the instability leads to the experimentally observed for- 
mation of asters and vortices. Similar vortices were ob- 
served in a system of interacting granular rods |26l l3l| . 
While we find that it suffices to consider only the density 
of tubules to explain the basic phenomenology, a better 
agreement with experiment is obtained when we include 
variable motor density and the motors attachment to the 
substrate HEEl. 



gleet the derivatives of m where it is appropriate, and 
the resulting equations assume the form: 



(75) 



(76) 
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FIG. 13: Motor contrectation profiles m for different values 
of (o for isolated aster for H = 0.125, B = 0.06, D = 5, and 
m = l,po = 4. 



Many aspects of self-assembly in tubule-motor systems 
require further inverstigation. In particular, we antic- 
ipate that flexibility of the microtubules may have a 
strong effect on the details of interaction, the question 
which we plan to address in future work. Another inter- 
esting question is role of the hydrodynamic interaction 
between the microtubules and effects of fluctuations on 
the orientation transition. Furthermore, our theory is 
derived in the limit of small density of the microtubules 
and takes into account only binary interactions among 
tubules. Certainly, multiparticle interactions are impor- 
tant in the high density state, such as bundles. General- 
ization of our work to the multiparticle interaction is a 
very challenging problem. 

There are many predictions following from our analy- 
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example, we find that the anisotropy of the interaction 
kernel is associated with the inhomogeneous density of 
motors along the microtubules. We observed that the 
motors attached to the substrate reduce the density of 
microtubules in the cores of asters and vortices. We pre- 
dicted an acceleration instability which leads to a drift 
of isolated asters. Hopefully, new generations of experi- 
ments will be able to address these issues. 

We thank Jacques Prost, Anthony Maggs, Leo 
Kadanoff and Valerii Vinokur for useful discussions. 
This work was supported by the U.S. Departemtn of 
FIG. 14: Color-coded images of density of motors at t = Ene ^ ^ s W-31-109-ENG-38 (IA) and DE-FG02- 
200 (left) and t = 1000 (right), red corresponds to maximum 04ER46135 (LT). 
density and blue to minimum density, domain size is 80 x 80 
units, C,o = 0.5, H — 0.125, D = 5, mo = 1. Arrows show 
orientation of microtubules. 




APPENDIX A: CONTINUUM EQUATIONS FOR 
ARBITRARY 7, O 



sis which possibly deserve experimental verification. For 



J 



We assume that two rods interact when angle between them is less or equal <fio independent of their spatial location. 
The master equation reads 



dPfat) _ D d 2 P(^t) 



dt 



+9 
+9 



i /2 r 2ir 



dtp 2 

2dydxP{x + y, t)P(x - y, t) [5(<j> - 7(3 - y) - (1 - -y)(x + y)) -5{4>-x- y)] 



n-4> /2 JO 



2tt 



2dydxP{x + y, t)P(x - y, t) [6(cf> - j(x - y) - (1 - <y)(x + y + 2tt)) - , 



x-y)] (Al) 



where x — (<pi + 4>2)/2 1 y — (<f>2 — 4>i)/2 (take into account that difiidfo — 2dxdy), and g is the collision rate, and we 
also included thermal diffusion of rods orientation oc D r . 
Performing integration over x from to 2-7T, we get 



dP{j>,t) d 2 P(^t) 

— D r — -tt-t, h 2g 



dt 



<Po/2 



dy [P(d> + 2 1V , t)P(<f> + 2(7 - l)y, t) - P{cj>, t)P(<t> - 2y, t)] 



+2g 



-0o/2 



dy [P{(j> + 2 ly + 2^7, t)P{cj> + 2(7 - l)y + 27T 7 , t) - P(cf>, t)P{cf> - 2y, t)} 



(A2) 



Changing variable y — > y + n in the second term in l|A2(l , we obtain 



g*W) n d 2 P{^t) 



4>o/2 



-</>o/2 



dy [P(<j> + 2y 7 , t)P(cj> + 2y( 7 - 1), t) - P(<j>, t)P(</> - 2y, t)} 



(A3) 



The case when all rods interact corresponds to (po = ir, and Ea . (|A3(I simplifies to 



dP{^t) d 2 P(^t) 
—^—+ gP{<j „t)=D r ——+2 9 



ir/2 



dy[P{<t> + 2y 1 ,t)P{(t> + 2y{ 1 -l),t) 



-tt/2 



(A4) 



By substitution y — > w/2 Eq. (| A3|) can be transformed into the form 
dP((/),t) d 2 P((j>,t) 



dt 



dw [P{<t> + vrf, t)P{<t> + w(r/ - 1), t) - P{<i>, t)P(4> - w, t)] 



(A5) 



where we changed variables t — > D r t, P — > gP/D r . 
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Let us consider a Fourier expansion of the probability distribution 

oo 

P(0,f) = p k(t)e ik * (A6) 



fe=-c 



where P_fc = P£. The Fourier harmonics P& are given by angular averages of exp(z/c0), see Eq. ©. The constant 
zeroth harmonic Pq = l/2itp, where p is the number density, 



P = / #P(0,r) =2nP (A7) 
Jo 

and the real and imaginary parts of Pi represent the components of the orientation vector r = ((cos0), (sine/))). 
Accordingly, t x + ir y = (exp[z0]) = Pf. 

After substitution of I|A6J| into Ea. (jA5jl we obtain the infinite series of equations for Pk 

Pk + k 2 P k = 20 o ^ PnP m (S[Mn-Y + Hi ~ 1))] " S{m^)]5 n+m , k (A8) 

n m 

where S{x) = sinx/x, and 5 n + m .k is the Kroneker symbol. 
For 0o = ft", the latter equation simplifies to 

P k + (fc 2 + l)P k = 2nJ2Y, PnPmS[n(n-f + m( 7 - l)]5 n+m>k (A9) 

n m 

Now we have to truncate this series. Assuming P n — for all \n\ > 2, from Ea. (|A8|l one obtains Pq = 0, 

A + Pi = P O Pi20o [S[Ml - 1)] + S{4>ol} - S(M - 1] 
+ 20 O P 2 P* [S[Ml + 1)] + S[<M7 - 2)] - S(2M - S(M (A10) 

and 

P 2 + 4P 2 = P O P 2 20 O [S[2Ml - 1)] + S[20 o7 ] - S(2<j> ) - 1] + 20 O P 2 [S[0 O (2 7 ~ 1)] " S(<M] (All) 
Neglecting the time derivative P2, we obtain P2 = j4P 2 with 

A= S[0o(2 7 - 1)] - S(M 

2/0o - (S[20 o (7 - 1)] + W07] - S(20 o ) - l)p/2n 1 j 

That allows us to close the equation for Pi, 

p + p = pp^ovr- 1 [s[Mi - 1)] + s[<M - s(M - 1] 

+2A0 o |P| 2 Pi [S[Ml + 1)] + S[Ml - 2)] - S(20o) - S(M] (A13) 

For 0o = 7T, 

P + (1 + p)P = pPa [5[tt(7 - 1)] + Shi] + 2 7 rA |P 1 | 2 P 1 [5[vr( 7 + 1)] + S[tt( 7 - 2)]] (AH) 

where 

A gTTgKgT-l)] 

4-0S[27r( 7 -l)]+S[27r 7 ]-l)p 1 J 



As seen from this equation, for < 7 < 1 we obtain an ordering instability which for large enough p > p c leads to a 
spontaneous alignment of filaments. 
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APPENDIX B: EVALUATION OF TERMS Z ,Z X 
1. Isotropic term Zq. 

We introduce new variable £ = r ri , and obtain after simple algebra 

Z a = I di f ° dwW {\£\)[P{r- £/2,<f> + w 7 , t)P{r + £/2,<t> + w(j-l),t) -P(r, <f,,t)P(r - £, <f> - w,t)](Bl) 

Now we assume that the probability distributions are smooth functions of spatial coordinates on the scale of the 
rod length I, and expand them near r, 

P(r + £, 0, t) = P(r, 0, t) + (€ • V)P(r, 0, t) + -(£ ■ V) 2 P(r, 0, t) + 0(£ 3 ), (B2) 

Performing integration over £ using kernel l|45l) . we get 

/<t>0 
dw [P(r, + w 7 , i)P(r, + ib(7 - 1), t) - P(r, 0, i)P(r, - t)] 

b 2 Z 100 

+ — / dw[-2VP(r,0 + w 7 ,t)VP(r,0 + w( 7 - l),i) (B3) 

16 J-<Po 

+V 2 P(r, + wj, t)P(r, (j> + w(j - 1), i) + P(r, + tu 7 , i)V 2 P(r, </> + u>( 7 - 1), i) -4P(r, 0, i)V 2 P(r, <f>-w, tj] 
Now we expand P's in the Fourier series over <fi, P(r,cf>,t) = ^2 n P n e ln ^ . The kth Fourier component of Z reads 
Zq = 2<t> \] 6 n+m! kPnPm[S[4>o(n'y + m( 7 - 1))] - S{m(j) )} 

n m 

+ ^ W,fe[(-2VP n VP m + V 2 P„P m + P„V 2 P m )S[n 7 o + m( 7 - 1)0 O ] - 4P„V 2 P m S(™0 o )] (B4) 



The first sum in Ea. (|B4|) coincides with the spatially uniform case (|A8|) . For k — we obtain (keeping only terms up 
to \n\, \m\ = 1) 

Z° = -^[V 2 (P 2 ) + 2V 2 (|P 1 | 2 )S(0 O )] (B5) 

For k = 1 

Zq 1 = 20 o PoP - 1)] + 5[0 o7 ] - S(<f> ) - 1] + 2^^* [Sfeofr + 1)] + Sfoofr - 2)] - S(20 o ) - Sfa)] 



8 



-[(V 2 P P! + P V 2 P - 2VPoVP!)(5[ 7 o ] + 5[( 7 - !)&,]) - 4P O V 2 P!S[0 O ] - 4P!V 2 P ] 



j 2 i 

+-P-[(-2VP?VP 2 + V 2 P*P 2 + P 1 *V 2 P 2 )5[(2 - 7 )0 O ] + (-2VP 2 VP* + V 2 P 2 Pi + P 2 V 2 P*)S[(1 + l)M 

o 

-4P*V 2 P 2 S(20 O ) - 4P 2 V 2 P 1 *S(0 O )] (B6) 

[this equations derived with the aid of Mathematica 32]]. 
We can again use the expression P 2 = AP 2 with constant A obtained for the spatially uniform case l|A12l) . 
Now, if we set 4>o — TTj 7 = 1/2 an d neglect higher order terms in the differential operators, we obtain 

Z° a = -^V 2 (P 2 ) 

Z\ = 2P Pi(4 - tt) - -P-P^ + y [V 2 (P Pi) - 4VP VP - ^PiV 2 P - {V 2 {P 2 Pt) - 4VP 1 *VP 2 )/3] 
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2. Anisotropic term Z\. 

For compactness of notations we introduce the following definitions z = x + iy,ip = 

= exp[i(j)]. Scalar 

product assumes the form r • n = Re(z* exp[i(f>]). Now Eq. I|58|) can be written in the form 



Zi 



dv ± j dw 



2W (2(n - r))Re(2( Zl - z )* {e^^ - e lW,+(7 - 1)to) )P(ri, </> + jw, t)P(2r - n, $ + (7 - l)w, t) 
-Wo(n - r)Re(( Zl - z)*(e^-^ - e«))P(i, 0, f)P(n, <j> - w, t)] (B7) 

Let us introduce £ = 2(ri — r), £ = 2(zi — z) in the first integrand of Ea. (|B7(l and £ = ri — r, £ — z\ — z in the 
second, and obtain 

r r<fro 
Z x = d£ dwWo(\m 

J J — <f) 

Pe(C*(e^ +7to) - e l (*+(T-i)»))p( r + £/ 2 , ^ + 7U ,, t)p( r - £/2, + (7 - l)w, t) 
- i?e(C*(e tW -"' ) - e**))P(r, 0, i)P(r + £, - to, t)l (B8) 

Now we again make use of the Fourier expansion P(r, <j>, t) = ^Pn^fje" 1 *. Multiplying Ea. HB8(l by (27r) _1 e~ jfc ^ 
and integrating over 4> from to 2-7T , we get 



2tt 



d«;Wb(|€l)e" 



Pe(C*(e lW+7U,) - e^+fr-^jPfr + £/2, + 7*0, t)P(r - £/2, + (7 - l)t«, i) 
- Pe(C*(e <( ^ w) - e l *))P(r, 0, i)P(r + £, - w, t) 



(B9) 



Using periodicity in </> we shift variables 



1 



jw in the first integrand: 



2tt 



dwW (\$\)e- ik *Re((*( 



[P(r + e/2, 0, t)P(r - |/2, - t)e ife ^ + P(r, 0, i)P(r + £, </> - w, t) 
Changing variables <f> — > <fi + w and then w — > —w where appropriate, we obtain 



(BIO) 



27T p4> 

d£ j d<t> dwW (\Z\)e- ik 't'Re((*e i 4') 



2ttj Jo 

P(r + g/2, 0, t)P(r - £/2, <f>-w, t)e ik ^ w - P(r + £/2, - to, t)P(r - £/2, 0, i)]e ife(1 -^ w 
+P(r, 0, i)P(r + £, - t) - P(r, - t)P(r + £, 0, i)e lto ] 
Let us consider again the zeroth and first moments only. For Z\ we get 

-1 p p2n p(f>o 

Z? = 27 J * / # y rf ^o(|l|)Pe(Ce^) 

[P(r + 1/2, 0, i)P(r - £/2, - w, t) - P(r + £/2, - tu, t)P(r - £/2, 0, i)] 
+P(r, 0, <)P(r + I, - w, t) - P(r, - «;, t)P(r + £, 0, t)] 

After cumbersome transformations in Mathcmatica it can be reduced to 

1 



(Bll) 



(B12) 



70 



64tt 7 



2tt /<<f>o 

d<f> / dw [3Re (V* exp(i^)) (P((/>)AP(<f> - w) - P{<f> - w)AP(4>)) 



21m (V* exp(i0)) (P(0) y P(0 - - P{<j)) x P{<l> ~ w) y )] 



(B13) 
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(here V = d x + id y , V* = d x — id y ). Performing integration over <fi and w, we obtain 



Z° = — 
1 32 



3V* (P_iAP - F AP-i) - - (V* (dyP-APo - d x P-id y P )) 



(B14) 



This expression can be written in the vector form 



1 16 







d 



3V • (rAp - pAr) - 2 — (d y T y d x p - d x T y d y p) + 2— (d y T x d x p - d x r x d y p) 



Similarly, after some algebra we also obtain the anisotropic part of the first moment, which for 7=1/2 reads 

Z\ = / dcj) dwe' 1 ^ [[Pe(e^V*)P(r, 0, t)P(r, - to, i) - P(r, 0, t)Pe(e i<#, V*)P(r, - to, i)]e"°/ 2 
4tt7 J_^ o L 

+P(r, 0, i)Pe(e^V*)P(r, - t) - P(r, - to, t)Pe(e^V*)P(r, ^, t)e™] (B15) 
Eq. (|B15|) can be written as 

b 2 ^ 



Z\ 



4tt 



hi J-<j>o 
ii»/2T> e (J<t>\ 



(B16) 



e™^Pe(e 4<p V*)P(r, 0, i)P(r, 0- to,t) + P(r, 0, i)Pe(e l0 V*)P(r, - to, i) 
Now performing integration over 0, to we obtain 

Z i = \~ t( 5 ((V2 - m)fo) - - m)0 o )) (V*P_ m P m + VP 2 _ m P ro ) 

m— 

+ (S(ro^o) - S((l/2 - m)&)) (V*P m P_ m + VP m P 2 _ ro )] 
Substituting </> = 7r, and keeping only the first two moments, we obtain 

Z l = [V*P P - Pi (V*P-1 + VP) + S 1/2 (V*P-lPl - V*PiP_i) + S 3/2 (V*i\P_l - V*P_lPl)] (B18) 

where S 1/2 = S{tt/2) and S 3/2 = S{3n/2). 
After some transformations we obtain 

■Kb 2 

Z\ = — [V*P P - (1 - Si/ 2 + S3/2) (V*P-iPi + VP1P1) - (S 1/2 - S3/2) (VP1P1 + V*PiP-i)] (B19) 



(B17) 



In the vector form it gives 



irb 2 



j^P^P ~ 2(1 - S 1/2 + S 3/2 )tV ■ t - 2(S 1/2 - S 3/2 )(tV)t 



Since S\/ 2 = 2/ir, S3/2 = — 2/37T, we obtain 



Z\ = b 2 



-pVp- T(V . T )--( T V)r 



(B20) 
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